#Hawkish Biases Replication 1
#Last modified: November 24, 2021
#Josh Kertzer

## This R file contains the code necessary to replicate the analysis in the main text. Its companion R file ("Hawkish Biases Replication 2.R") replicates the analysis in the supplementary appendix.
# All of the following analyses were carried out using R version 4.1.2. GUI 1.77 Big Sur Arm Build (8007) on an M1 Pro Macbook Pro running MacOS Monterey; they were also replicated on a 2017 3 Ghz 10-Core Intel Xeon W iMac Pro. 

#Load libraries
library(here)
library(dplyr)
library(ggplot2)
library(stargazer)
library(sparsereg)
library(caret)
library(caretEnsemble)
library(tidyr)

#Load data
dat <- get(load("Hawkish-Biases.RData"))

#Load functions

#Calculate mode
calc.mode <- function(x){
	unique.x <- unique(x)
	unique.x[which.max(tabulate(match(x, unique.x)))]
}

#Rescale variables from 0-1
rescale <- function(x){
	return((x-min(x,na.rm=TRUE))/(max(x-min(x,na.rm=TRUE),na.rm=TRUE)))
}

#Return confidence intervals
ateCalc <- function(mod, x=2){
	b <- coef(mod)[x]
	s <- sqrt(vcov(mod)[x,x])
	return(c(b,b-1.96*s,b+1.96*s, b-1.65*s,b+1.65*s)) 
}


#### Build datasets for analysis
dat$groupCond <- as.character(dat$groupCond)
dat.i <- dat[which(dat$groupCond=="Individual"),] #Individual condition
dat.f <- dat[which(dat$groupCond=="Horizontal"),] #Horizontal condition
dat.h <- dat[which(dat$groupCond=="Hierarchical"),] #Hierarchical condition

#Create different versions of group conditions

### Hierarchical ##
dat.h.1 <- dat[which(dat$groupCond=="Hierarchical" & dat$role_pt=="Leader"),] #Hierarchical condition

#But, also create dissensus measure, capture traits, for hierarchical groups
adv.h <- dat.h[which(dat.h$role_pt!="Leader"),] #Advisor-specific dataframe

#Diversity dataframes:
df_dem <- cbind(rescale(dat.h$age), dat.h$male, dat.h$educ1, dat.h$income1, dat.h$christian, dat.h$white) #Demographic
df_att <- cbind(dat.h$sdo1, dat.h$rwa1, dat.h$mi1, dat.h$iso1, dat.h$ideo1) #Political attitudes
df_disp <- cbind(dat.h$cognition1, dat.h$persExtra1, dat.h$persAgree1, dat.h$persConsc1, dat.h$persNeuro1, dat.h$persOpen1, dat.h$aggression1, dat.h$risk1) #Dispositional/personality charactistics
df_exp <- cbind(rescale(log(dat.h$maxemp+0.01)), rescale(dat.h$yrsworked), dat.h$smallgrp1, dat.h$govemp1) #Experience

## Prospect theory experiment
dat.h.1.pt <- dat.h.1
dat.h.1.pt$dissensus <- NA #Dissensus
dat.h.1.pt$diversity_dem <- dat.h.1.pt$diversity_att <- dat.h.1.pt$diversity_exp <- dat.h.1.pt$diversity_disp <- NA #Diversity within hierarchical group as a whole

for (i in unique(dat.h.1$GroupID_pt)){
	j <- which(dat.h$GroupID_pt==i) #Full dataframe
	k <- which(dat.h.1.pt$GroupID_pt==i) #One observation per group (for leaders)
	l <- which(adv.h$GroupID_pt==i) #Advisors only
	
	#Diversity calculations
	dat.h.1.pt$diversity_dem[k] <- mean(apply(df_dem[j,],2,var,na.rm=TRUE))
	dat.h.1.pt$diversity_att[k] <- mean(apply(df_att[j,],2,var,na.rm=TRUE))
	dat.h.1.pt$diversity_disp[k] <- mean(apply(df_disp[j,],2,var,na.rm=TRUE))
	dat.h.1.pt$diversity_exp[k] <- mean(apply(df_exp[j,],2,var,na.rm=TRUE))
	#Dissensus across the group as a whole
	dat.h.1.pt$dissensus[k] <- var(dat.h$riskyDV_pt[j], na.rm=TRUE)
	#Demographics  - utilize whole group composition
	dat.h.1.pt$mi1_g[k] <- mean(dat.h$mi1[j],na.rm=TRUE)
	dat.h.1.pt$iso1_g[k] <- mean(dat.h$iso1[j],na.rm=TRUE)
	dat.h.1.pt$risk1_g[k] <- mean(dat.h$risk1[j],na.rm=TRUE)
	dat.h.1.pt$sdo1_g[k] <- mean(dat.h$sdo1[j],na.rm=TRUE)
	dat.h.1.pt$rwa1_g[k] <- mean(dat.h$rwa1[j],na.rm=TRUE)
	dat.h.1.pt$persExtra1_g[k] <- mean(dat.h$persExtra1[j],na.rm=TRUE)		
	dat.h.1.pt$persAgree1_g[k] <- mean(dat.h$persAgree1[j],na.rm=TRUE)
	dat.h.1.pt$persConsc1_g[k] <- mean(dat.h$persConsc1[j],na.rm=TRUE)
	dat.h.1.pt$persNeuro1_g[k] <- mean(dat.h$persNeuro1[j],na.rm=TRUE)
	dat.h.1.pt$persOpen1_g[k] <- mean(dat.h$persOpen1[j],na.rm=TRUE)
	dat.h.1.pt$cognition1_g[k] <- mean(dat.h$cognition1[j],na.rm=TRUE)
	dat.h.1.pt$aggression1_g[k] <- mean(dat.h$aggression1[j],na.rm=TRUE)
	dat.h.1.pt$age_g[k] <- mean(dat.h$age[j],na.rm=TRUE)
	dat.h.1.pt$male_g[k] <- mean(dat.h$male[j],na.rm=TRUE)	
	dat.h.1.pt$christian_g[k] <- mean(dat.h$christian[j],na.rm=TRUE)
	dat.h.1.pt$white_g[k] <- mean(dat.h$white[j],na.rm=TRUE)
	dat.h.1.pt$pid1_g[k] <- mean(dat.h$pid1[j],na.rm=TRUE)
	dat.h.1.pt$ideo1_g[k] <- mean(dat.h$ideo1[j],na.rm=TRUE)
	dat.h.1.pt$interest1_g[k] <- mean(dat.h$interest1[j],na.rm=TRUE)
	dat.h.1.pt$educ1_g[k] <- mean(dat.h$educ1[j],na.rm=TRUE)
	dat.h.1.pt$income1_g[k] <- mean(dat.h$income1[j],na.rm=TRUE)
	#Level of group participation (per capita)
	dat.h.1.pt$timesSpoken_pt_g[k] <- sum(dat.h$timesSpoken_pt[j], na.rm=TRUE)
	dat.h.1.pt$timesSpokenPC_pt_g[k] <- dat.h.1.pt$timesSpoken_pt_g[k]/dat.h.1.pt$groupN_pt[k]
}


## Intentionality bias experiment
dat.h.1.ib <- dat.h.1
dat.h.1.ib$dissensus <- NA #Dissensus
dat.h.1.ib$diversity_dem <- dat.h.1.ib$diversity_att <- dat.h.1.ib$diversity_exp <- dat.h.1.ib$diversity_disp <- NA #Diversity within hierarchical group as a whole

for (i in unique(dat.h.1$GroupID_ib)){
	j <- which(dat.h$GroupID_ib==i) #Full dataframe
	k <- which(dat.h.1.ib$GroupID_ib==i) #One observation per group (for leaders)
	l <- which(adv.h$GroupID_ib==i) #Advisors only
	
	#Diversity calculations
	dat.h.1.ib$diversity_dem[k] <- mean(apply(df_dem[j,],2,var,na.rm=TRUE))
	dat.h.1.ib$diversity_att[k] <- mean(apply(df_att[j,],2,var,na.rm=TRUE))
	dat.h.1.ib$diversity_disp[k] <- mean(apply(df_disp[j,],2,var,na.rm=TRUE))
	dat.h.1.ib$diversity_exp[k] <- mean(apply(df_exp[j,],2,var,na.rm=TRUE))	
	#Dissensus across the group as a whole
	dat.h.1.ib$dissensus[k] <- var(dat.h$intentDV_ib[j], na.rm=TRUE)
	#Demographics  - assume whole group composition (already have same variables saved as individuals, without the _g)
	dat.h.1.ib$mi1_g[k] <- mean(dat.h$mi1[j],na.rm=TRUE)
	dat.h.1.ib$iso1_g[k] <- mean(dat.h$iso1[j],na.rm=TRUE)
	dat.h.1.ib$risk1_g[k] <- mean(dat.h$risk1[j],na.rm=TRUE)
	dat.h.1.ib$sdo1_g[k] <- mean(dat.h$sdo1[j],na.rm=TRUE)
	dat.h.1.ib$rwa1_g[k] <- mean(dat.h$rwa1[j],na.rm=TRUE)
	dat.h.1.ib$persExtra1_g[k] <- mean(dat.h$persExtra1[j],na.rm=TRUE)		
	dat.h.1.ib$persAgree1_g[k] <- mean(dat.h$persAgree1[j],na.rm=TRUE)
	dat.h.1.ib$persConsc1_g[k] <- mean(dat.h$persConsc1[j],na.rm=TRUE)
	dat.h.1.ib$persNeuro1_g[k] <- mean(dat.h$persNeuro1[j],na.rm=TRUE)
	dat.h.1.ib$persOpen1_g[k] <- mean(dat.h$persOpen1[j],na.rm=TRUE)
	dat.h.1.ib$cognition1_g[k] <- mean(dat.h$cognition1[j],na.rm=TRUE)
	dat.h.1.ib$aggression1_g[k] <- mean(dat.h$aggression1[j],na.rm=TRUE)
	dat.h.1.ib$age_g[k] <- mean(dat.h$age[j],na.rm=TRUE)
	dat.h.1.ib$male_g[k] <- mean(dat.h$male[j],na.rm=TRUE)	
	dat.h.1.ib$christian_g[k] <- mean(dat.h$christian[j],na.rm=TRUE)
	dat.h.1.ib$white_g[k] <- mean(dat.h$white[j],na.rm=TRUE)
	dat.h.1.ib$pid1_g[k] <- mean(dat.h$pid1[j],na.rm=TRUE)
	dat.h.1.ib$ideo1_g[k] <- mean(dat.h$ideo1[j],na.rm=TRUE)
	dat.h.1.ib$interest1_g[k] <- mean(dat.h$interest1[j],na.rm=TRUE)
	dat.h.1.ib$educ1_g[k] <- mean(dat.h$educ1[j],na.rm=TRUE)
	dat.h.1.ib$income1_g[k] <- mean(dat.h$income1[j],na.rm=TRUE)	
	#Level of group participation (per capita)
	dat.h.1.ib$timesSpoken_ib_g[k] <- sum(dat.h$timesSpoken_ib[j], na.rm=TRUE)
	dat.h.1.ib$timesSpokenPC_ib_g[k] <- dat.h.1.ib$timesSpoken_ib_g[k]/dat.h.1.ib$groupN_ib[k]	
}

## Reactive devaluation experiment

#RD submodule
dat.h.1.rd <- dat.h.1
dat.h.1.rd$dissensus <- NA #Dissensus
dat.h.1.rd$diversity_dem <- dat.h.1.rd$diversity_att <- dat.h.1.rd$diversity_exp  <- dat.h.1.rd$diversity_disp <- NA #Diversity within hierarchical group as a whole

for (i in unique(dat.h.1$GroupID_rd)){
	j <- which(dat.h$GroupID_rd==i) #Full dataframe
	k <- which(dat.h.1.rd$GroupID_rd==i) #One observation per group (for leaders)
	l <- which(adv.h$GroupID_rd==i) #Advisors only
	
	#Diversity calculations
	dat.h.1.rd$diversity_dem[k] <- mean(apply(df_dem[j,],2,var,na.rm=TRUE))
	dat.h.1.rd$diversity_att[k] <- mean(apply(df_att[j,],2,var,na.rm=TRUE))
	dat.h.1.rd$diversity_exp[k] <- mean(apply(df_exp[j,],2,var,na.rm=TRUE))
	dat.h.1.rd$diversity_disp[k] <- mean(apply(df_disp[j,],2,var,na.rm=TRUE))
	#Dissensus across the group as a whole
	dat.h.1.rd$dissensus[k] <- var(dat.h$supportDV_rd[j], na.rm=TRUE)
	#Demographics  - assume whole group composition (already have some variables saved as individuals, without the _g)
	dat.h.1.rd$mi1_g[k] <- mean(dat.h$mi1[j],na.rm=TRUE)
	dat.h.1.rd$iso1_g[k] <- mean(dat.h$iso1[j],na.rm=TRUE)
	dat.h.1.rd$risk1_g[k] <- mean(dat.h$risk1[j],na.rm=TRUE)
	dat.h.1.rd$sdo1_g[k] <- mean(dat.h$sdo1[j],na.rm=TRUE)
	dat.h.1.rd$rwa1_g[k] <- mean(dat.h$rwa1[j],na.rm=TRUE)
	dat.h.1.rd$persExtra1_g[k] <- mean(dat.h$persExtra1[j],na.rm=TRUE)		
	dat.h.1.rd$persAgree1_g[k] <- mean(dat.h$persAgree1[j],na.rm=TRUE)
	dat.h.1.rd$persConsc1_g[k] <- mean(dat.h$persConsc1[j],na.rm=TRUE)
	dat.h.1.rd$persNeuro1_g[k] <- mean(dat.h$persNeuro1[j],na.rm=TRUE)
	dat.h.1.rd$persOpen1_g[k] <- mean(dat.h$persOpen1[j],na.rm=TRUE)
	dat.h.1.rd$cognition1_g[k] <- mean(dat.h$cognition1[j],na.rm=TRUE)
	dat.h.1.rd$aggression1_g[k] <- mean(dat.h$aggression1[j],na.rm=TRUE)
	dat.h.1.rd$age_g[k] <- mean(dat.h$age[j],na.rm=TRUE)
	dat.h.1.rd$male_g[k] <- mean(dat.h$male[j],na.rm=TRUE)	
	dat.h.1.rd$christian_g[k] <- mean(dat.h$christian[j],na.rm=TRUE)
	dat.h.1.rd$white_g[k] <- mean(dat.h$white[j],na.rm=TRUE)
	dat.h.1.rd$pid1_g[k] <- mean(dat.h$pid1[j],na.rm=TRUE)
	dat.h.1.rd$ideo1_g[k] <- mean(dat.h$ideo1[j],na.rm=TRUE)
	dat.h.1.rd$interest1_g[k] <- mean(dat.h$interest1[j],na.rm=TRUE)
	dat.h.1.rd$educ1_g[k] <- mean(dat.h$educ1[j],na.rm=TRUE)
	dat.h.1.rd$income1_g[k] <- mean(dat.h$income1[j],na.rm=TRUE)			
	#Level of group participation (per capita)
	dat.h.1.rd$timesSpoken_rd_g[k] <- sum(dat.h$timesSpoken_rd[j], na.rm=TRUE)
	dat.h.1.rd$timesSpokenPC_rd_g[k] <- dat.h.1.rd$timesSpoken_rd_g[k]/dat.h.1.rd$groupN_rd[k]	
}


### Horizontal conditions 

#Create three decision rules: 
#Unanimity: all group members generate same answer
#Majority rule: a majority of respondents in a gorup of a given size give the same answer
#Median voter: the decision rule used for the main analysis

#Diversity dataframes (use for all submodules):
df_dem <- cbind(rescale(dat.f$age), dat.f$male, dat.f$educ1, dat.f$income1, dat.f$christian, dat.f$white) #Demographic
df_att <- cbind(dat.f$sdo1, dat.f$rwa1, dat.f$mi1, dat.f$iso1, dat.f$ideo1) #Political attitudes
df_disp <- cbind(dat.f$cognition1, dat.f$persExtra1, dat.f$persAgree1, dat.f$persConsc1, dat.f$persNeuro1, dat.f$persOpen1, dat.f$aggression1, dat.f$risk1) #Dispositional/personality charactistics
df_exp <- cbind(rescale(log(dat.f$maxemp+0.01)), rescale(dat.f$yrsworked), dat.f$smallgrp1, dat.f$govemp1) #Experience

#Note: we need to create horizontal dataframes separately using submodule-level group IDs

## Prospect theory experiment
dat.f.1.pt <- data.frame(GroupID_pt=unique(dat.f$GroupID_pt), groupCond="Horizontal", stringsAsFactors=FALSE)

dat.f.1.pt$un_grp_pt <- NA #Strict unanimity
dat.f.1.pt$maj_grp_pt <-  NA #Majority rule
dat.f.1.pt$riskyDV_pt <- NA #Group decision (mode)
dat.f.1.pt$medDV_pt <- NA #Median vote 
dat.f.1.pt$gainTrt_pt <- NA #Treatment
dat.f.1.pt$lossTrt_pt <- NA #Treatment
dat.f.1.pt$groupN_pt <- NA #Group size
dat.f.1.pt$risk_m_pt <- dat.f.1.pt$risk_mw_pt <- dat.f.1.pt$risk_min_pt <- dat.f.1.pt$risk_max_pt <- NA #Risk attitudes
dat.f.1.pt$diversity_dem <- dat.f.1.pt$diversity_att <- dat.f.1.pt$diversity_exp <- NA #Group diversity measures (demographic, attitudinal, experience)
dat.f.1.pt$dissensus <- NA #Dissensus measure (variance of DV in group)
dat.f.1.pt$numF <- dat.f.1.pt$pctF <- NA #Number of women in group, and % of women in group

for (i in unique(dat.f$GroupID_pt)){
	j <- which(dat.f$GroupID_pt==i)
	k <- which(dat.f.1.pt$GroupID_pt==i)
	
	modal.answer <- calc.mode(dat.f$riskyDV_pt[j])
	num.mode <- length(which(dat.f$riskyDV_pt[j]==modal.answer))
	group.size <- length(na.omit(dat.f$riskyDV_pt[j]))
	dat.f.1.pt$groupN_pt[k] <- group.size	
	if (group.size >=3){
		dat.f.1.pt$un_grp_pt[k] <- ifelse(num.mode==group.size,1,0)
		dat.f.1.pt$maj_grp_pt[k] <- ifelse(num.mode>(group.size/2),1,0)
		dat.f.1.pt$riskyDV_pt[k] <- modal.answer
		dat.f.1.pt$medDV_pt[k] <- median(dat.f$riskyDV_pt[j], na.rm=TRUE)
		dat.f.1.pt$gainTrt_pt[k] <- ifelse(length(unique(dat.f$gainTrt_pt[j]))==1,unique(dat.f$gainTrt_pt[j]),NA)
		dat.f.1.pt$lossTrt_pt[k] <- ifelse(length(unique(dat.f$lossTrt_pt[j]))==1,unique(dat.f$lossTrt_pt[j]),NA)		
		dat.f.1.pt$risk_m_pt[k] <- mean(dat.f$risk1[j],na.rm=TRUE)
		dat.f.1.pt$risk_mw_pt[k] <- mean(dat.f$risk1[j]*log(dat.f$wordsSpoken_pt[j]+0.01),na.rm=TRUE)		
		dat.f.1.pt$risk_min_pt[k] <- min(dat.f$risk1[j],na.rm=TRUE)		
		dat.f.1.pt$risk_max_pt[k] <- max(dat.f$risk1[j],na.rm=TRUE)	
		#Diversity calculations
		dat.f.1.pt$diversity_dem[k] <- mean(apply(df_dem[j,],2,var,na.rm=TRUE))
		dat.f.1.pt$diversity_att[k] <- mean(apply(df_att[j,],2,var,na.rm=TRUE))
		dat.f.1.pt$diversity_disp[k] <- mean(apply(df_disp[j,],2,var,na.rm=TRUE))
		dat.f.1.pt$diversity_exp[k] <- mean(apply(df_exp[j,],2,var,na.rm=TRUE))
		#Dissensus
		dat.f.1.pt$dissensus[k] <- var(dat.f$riskyDV_pt[j],na.rm=TRUE)
		#Demographics  - assume whole group composition
		dat.f.1.pt$mi1[k] <- mean(dat.f$mi1[j],na.rm=TRUE)
		dat.f.1.pt$iso1[k] <- mean(dat.f$iso1[j],na.rm=TRUE)
		dat.f.1.pt$risk1[k] <- mean(dat.f$risk1[j],na.rm=TRUE)
		dat.f.1.pt$sdo1[k] <- mean(dat.f$sdo1[j],na.rm=TRUE)
		dat.f.1.pt$rwa1[k] <- mean(dat.f$rwa1[j],na.rm=TRUE)
		dat.f.1.pt$persExtra1[k] <- mean(dat.f$persExtra1[j],na.rm=TRUE)		
		dat.f.1.pt$persAgree1[k] <- mean(dat.f$persAgree1[j],na.rm=TRUE)
		dat.f.1.pt$persConsc1[k] <- mean(dat.f$persConsc1[j],na.rm=TRUE)
		dat.f.1.pt$persNeuro1[k] <- mean(dat.f$persNeuro1[j],na.rm=TRUE)
		dat.f.1.pt$persOpen1[k] <- mean(dat.f$persOpen1[j],na.rm=TRUE)
		dat.f.1.pt$cognition1[k] <- mean(dat.f$cognition1[j],na.rm=TRUE)
		dat.f.1.pt$aggression1[k] <- mean(dat.f$aggression1[j],na.rm=TRUE)
		dat.f.1.pt$age[k] <- mean(dat.f$age[j],na.rm=TRUE)
		dat.f.1.pt$male[k] <- mean(dat.f$male[j],na.rm=TRUE)			
		dat.f.1.pt$christian[k] <- mean(dat.f$christian[j],na.rm=TRUE)
		dat.f.1.pt$white[k] <- mean(dat.f$white[j],na.rm=TRUE)
		dat.f.1.pt$pid1[k] <- mean(dat.f$pid1[j],na.rm=TRUE)
		dat.f.1.pt$ideo1[k] <- mean(dat.f$ideo1[j],na.rm=TRUE)
		dat.f.1.pt$interest1[k] <- mean(dat.f$interest1[j],na.rm=TRUE)
		dat.f.1.pt$educ1[k] <- mean(dat.f$educ1[j],na.rm=TRUE)
		dat.f.1.pt$income1[k] <- mean(dat.f$income1[j],na.rm=TRUE)	
		#Gender composition
		dat.f.1.pt$numF[k] <- length(which(dat.f$male[j]==0))
		dat.f.1.pt$pctF[k] <- (dat.f.1.pt$numF[k]/group.size)
		#Level of group participation (per capita)
		dat.f.1.pt$timesSpoken_pt_g[k] <- sum(dat.f$timesSpoken_pt[j], na.rm=TRUE)
		dat.f.1.pt$timesSpokenPC_pt_g[k] <- dat.f.1.pt$timesSpoken_pt_g[k]/dat.f.1.pt$groupN_pt[k]			
	}
}

## Intentionality bias experiment

dat.f.1.ib <- data.frame(GroupID_ib=unique(dat.f$GroupID_ib), groupCond="Horizontal", stringsAsFactors=FALSE)

dat.f.1.ib$un_grp_ib <- NA #Strict unanimity
dat.f.1.ib$maj_grp_ib <- NA #Majority rule 
dat.f.1.ib$intentDV_ib <- NA #Group decision (mode)
dat.f.1.ib$medDV_ib <- NA #Median voter
dat.f.1.ib$fatalTrt_ib <- NA #Treatments
dat.f.1.ib$groupN_ib <- NA #Group size
dat.f.1.ib$MI_m_ib <- dat.f.1.ib$MI_mw_ib <- dat.f.1.ib$MI_min_ib <- dat.f.1.ib$MI_max_ib <- NA #Hawkishness
dat.f.1.ib$agg_m_ib <- dat.f.1.ib$agg_mw_ib <- dat.f.1.ib$agg_min_ib <- dat.f.1.ib$agg_max_ib <- NA #Aggression
dat.f.1.ib$diversity_dem <- dat.f.1.ib$diversity_att <- dat.f.1.ib$diversity_exp <- dat.f.1.ib$diversity_disp <- NA #Group diversity measures (demographic, attitudinal, experience)
dat.f.1.ib$dissensus <- NA #Dissensus measure (variance of DV in group)
dat.f.1.ib$dissensus <- NA #Dissensus measure (variance of DV in group)
dat.f.1.ib$numF <- dat.f.1.ib$pctF <- NA #Number of women in group, and % of women in group


for (i in unique(dat.f$GroupID_ib)){
	j <- which(dat.f$GroupID_ib==i)
	k <- which(dat.f.1.ib$GroupID_ib==i)
	
	modal.answer <- calc.mode(dat.f$intentDV_ib[j])
	num.mode <- length(which(dat.f$intentDV_ib[j]==modal.answer))
	group.size <- length(na.omit(dat.f$intentDV_ib[j]))
	dat.f.1.ib$groupN_ib[k] <- group.size
	if (group.size >= 3){
		dat.f.1.ib$un_grp_ib[k] <- ifelse(num.mode==group.size,1,0)
		dat.f.1.ib$maj_grp_ib[k] <- ifelse(num.mode>(group.size/2),1,0)
		dat.f.1.ib$intentDV_ib[k] <- modal.answer
		dat.f.1.ib$medDV_ib[k] <- median(dat.f$intentDV_ib[j], na.rm=TRUE)		
		dat.f.1.ib$fatalTrt_ib[k] <- ifelse(length(unique(dat.f$fatalTrt_ib[j]))==1,unique(dat.f$fatalTrt_ib[j]),NA)
		dat.f.1.ib$MI_m_ib[k] <- mean(dat.f$mi1[j],na.rm=TRUE)
		dat.f.1.ib$MI_mw_ib[k] <- mean(dat.f$mi1[j]*log(dat.f$wordsSpoken_ib[j]+0.01),na.rm=TRUE)		
		dat.f.1.ib$MI_min_ib[k] <- min(dat.f$mi1[j],na.rm=TRUE)		
		dat.f.1.ib$MI_max_ib[k] <- max(dat.f$mi1[j],na.rm=TRUE)
		dat.f.1.ib$agg_m_ib[k] <- mean(dat.f$aggression1[j],na.rm=TRUE)
		dat.f.1.ib$agg_mw_ib[k] <- mean(dat.f$aggression1[j]*log(dat.f$wordsSpoken_ib[j]+0.01),na.rm=TRUE)		
		dat.f.1.ib$agg_min_ib[k] <- min(dat.f$aggression1[j],na.rm=TRUE)		
		dat.f.1.ib$agg_max_ib[k] <- max(dat.f$aggression1[j],na.rm=TRUE)
		#Diversity calculations
		dat.f.1.ib$diversity_dem[k] <- mean(apply(df_dem[j,],2,var,na.rm=TRUE))
		dat.f.1.ib$diversity_att[k] <- mean(apply(df_att[j,],2,var,na.rm=TRUE))
		dat.f.1.ib$diversity_disp[k] <- mean(apply(df_disp[j,],2,var,na.rm=TRUE))
		dat.f.1.ib$diversity_exp[k] <- mean(apply(df_exp[j,],2,var,na.rm=TRUE))
		#Dissensus
		dat.f.1.ib$dissensus[k] <- var(dat.f$intentDV_ib[j],na.rm=TRUE)
		#Demographics  - assume whole group composition
		dat.f.1.ib$mi1[k] <- mean(dat.f$mi1[j],na.rm=TRUE)
		dat.f.1.ib$iso1[k] <- mean(dat.f$iso1[j],na.rm=TRUE)
		dat.f.1.ib$risk1[k] <- mean(dat.f$risk1[j],na.rm=TRUE)
		dat.f.1.ib$sdo1[k] <- mean(dat.f$sdo1[j],na.rm=TRUE)
		dat.f.1.ib$rwa1[k] <- mean(dat.f$rwa1[j],na.rm=TRUE)
		dat.f.1.ib$persExtra1[k] <- mean(dat.f$persExtra1[j],na.rm=TRUE)		
		dat.f.1.ib$persAgree1[k] <- mean(dat.f$persAgree1[j],na.rm=TRUE)
		dat.f.1.ib$persConsc1[k] <- mean(dat.f$persConsc1[j],na.rm=TRUE)
		dat.f.1.ib$persNeuro1[k] <- mean(dat.f$persNeuro1[j],na.rm=TRUE)
		dat.f.1.ib$persOpen1[k] <- mean(dat.f$persOpen1[j],na.rm=TRUE)
		dat.f.1.ib$cognition1[k] <- mean(dat.f$cognition1[j],na.rm=TRUE)
		dat.f.1.ib$aggression1[k] <- mean(dat.f$aggression1[j],na.rm=TRUE)
		dat.f.1.ib$age[k] <- mean(dat.f$age[j],na.rm=TRUE)
		dat.f.1.ib$male[k] <- mean(dat.f$male[j],na.rm=TRUE)	
		dat.f.1.ib$christian[k] <- mean(dat.f$christian[j],na.rm=TRUE)
		dat.f.1.ib$white[k] <- mean(dat.f$white[j],na.rm=TRUE)
		dat.f.1.ib$pid1[k] <- mean(dat.f$pid1[j],na.rm=TRUE)
		dat.f.1.ib$ideo1[k] <- mean(dat.f$ideo1[j],na.rm=TRUE)
		dat.f.1.ib$interest1[k] <- mean(dat.f$interest1[j],na.rm=TRUE)
		dat.f.1.ib$educ1[k] <- mean(dat.f$educ1[j],na.rm=TRUE)
		dat.f.1.ib$income1[k] <- mean(dat.f$income1[j],na.rm=TRUE)	
		#Gender composition
		dat.f.1.ib$numF[k] <- length(which(dat.f$male[j]==0))
		dat.f.1.ib$pctF[k] <- (dat.f.1.ib$numF[k]/group.size)
		#Level of group participation (per capita)
		dat.f.1.ib$timesSpoken_ib_g[k] <- sum(dat.f$timesSpoken_ib[j], na.rm=TRUE)
		dat.f.1.ib$timesSpokenPC_ib_g[k] <- dat.f.1.ib$timesSpoken_ib_g[k]/dat.f.1.ib$groupN_ib[k]			
	}
}

## Reactive devaluation experiment
dat.f.1.rd <- data.frame(GroupID_rd=unique(dat.f$GroupID_rd), groupCond="Horizontal", stringsAsFactors=FALSE)

dat.f.1.rd$un_grp_rd <- NA #Strict unanimity
dat.f.1.rd$maj_grp_rd <- NA #Majority rule
dat.f.1.rd$supportDV_rd <- NA #Group decision (mode)
dat.f.1.rd$medDV_rd <- NA #Median voter
dat.f.1.rd$chinaTrt_rd <- NA #Treatments
dat.f.1.rd$groupN_rd <- NA #Group size
dat.f.1.rd$MI_m_rd <- dat.f.1.rd$MI_mw_rd <- dat.f.1.rd$MI_min_rd <- dat.f.1.rd$MI_max_rd <- NA #Hawkishness
dat.f.1.rd$agg_m_rd <- dat.f.1.rd$agg_mw_rd <- dat.f.1.rd$agg_min_rd <- dat.f.1.rd$agg_max_rd <- NA #Aggression
dat.f.1.rd$pid_m_rd <- dat.f.1.rd$pid_mw_rd <- dat.f.1.rd$pid_min_rd <- dat.f.1.rd$pid_max_rd <- NA #Party ID
dat.f.1.rd$diversity_dem <- dat.f.1.rd$diversity_att <- dat.f.1.rd$diversity_exp <-  dat.f.1.rd$diversity_disp <- NA #Group diversity measures (demographic, attitudinal, experience)
dat.f.1.rd$dissensus <- NA #Dissensus measure (variance of DV in group)
dat.f.1.rd$numF <- dat.f.1.rd$pctF <- NA #Number of women in group, and % of women in group


for (i in unique(dat.f$GroupID_rd)){
	j <- which(dat.f$GroupID_rd==i)
	k <- which(dat.f.1.rd$GroupID_rd==i)

	modal.answer <- calc.mode(dat.f$supportDV_rd[j])
	num.mode <- length(which(dat.f$supportDV_rd[j]==modal.answer))
	group.size <- length(na.omit(dat.f$supportDV_rd[j]))
	dat.f.1.rd$groupN_rd[k] <- group.size
	if (group.size >=3){
		dat.f.1.rd$un_grp_rd[k] <- ifelse(num.mode==group.size,1,0)
		dat.f.1.rd$maj_grp_rd[k] <- ifelse(num.mode>(group.size/2),1,0)
		dat.f.1.rd$supportDV_rd[k] <- modal.answer
		dat.f.1.rd$medDV_rd[k] <- median(dat.f$supportDV_rd[j],na.rm=TRUE)			
		dat.f.1.rd$chinaTrt_rd[k] <- ifelse(length(unique(dat.f$chinaTrt_rd[j]))==1,unique(dat.f$chinaTrt_rd[j]),NA)
		dat.f.1.rd$MI_m_rd[k] <- mean(dat.f$mi1[j],na.rm=TRUE)
		dat.f.1.rd$MI_mw_rd[k] <- mean(dat.f$mi1[j]*log(dat.f$wordsSpoken_rd[j]+0.01),na.rm=TRUE)		
		dat.f.1.rd$MI_min_rd[k] <- min(dat.f$mi1[j],na.rm=TRUE)		
		dat.f.1.rd$MI_max_rd[k] <- max(dat.f$mi1[j],na.rm=TRUE)
		dat.f.1.rd$agg_m_rd[k] <- mean(dat.f$aggression1[j],na.rm=TRUE)
		dat.f.1.rd$agg_mw_rd[k] <- mean(dat.f$aggression1[j]*log(dat.f$wordsSpoken_rd[j]+0.01),na.rm=TRUE)		
		dat.f.1.rd$agg_min_rd[k] <- min(dat.f$aggression1[j],na.rm=TRUE)		
		dat.f.1.rd$agg_max_rd[k] <- max(dat.f$aggression1[j],na.rm=TRUE)
		dat.f.1.rd$pid_m_rd[k] <- mean(dat.f$pid1[j],na.rm=TRUE)
		dat.f.1.rd$pid_mw_rd[k] <- mean(dat.f$pid1[j]*log(dat.f$wordsSpoken_rd[j]+0.01),na.rm=TRUE)		
		dat.f.1.rd$pid_min_rd[k] <- min(dat.f$pid1[j],na.rm=TRUE)		
		dat.f.1.rd$pid_max_rd[k] <- max(dat.f$pid1[j],na.rm=TRUE)
		#Diversity calculations
		dat.f.1.rd$diversity_dem[k] <- mean(apply(df_dem[j,],2,var,na.rm=TRUE))
		dat.f.1.rd$diversity_att[k] <- mean(apply(df_att[j,],2,var,na.rm=TRUE))
		dat.f.1.rd$diversity_disp[k] <- mean(apply(df_disp[j,],2,var,na.rm=TRUE))
		dat.f.1.rd$diversity_exp[k] <- mean(apply(df_exp[j,],2,var,na.rm=TRUE))
		#Dissensus
		dat.f.1.rd$dissensus[k] <- var(dat.f$supportDV_rd[j],na.rm=TRUE)
		#Demographics  - assume whole group composition
		dat.f.1.rd$mi1[k] <- mean(dat.f$mi1[j],na.rm=TRUE)
		dat.f.1.rd$iso1[k] <- mean(dat.f$iso1[j],na.rm=TRUE)
		dat.f.1.rd$risk1[k] <- mean(dat.f$risk1[j],na.rm=TRUE)
		dat.f.1.rd$sdo1[k] <- mean(dat.f$sdo1[j],na.rm=TRUE)
		dat.f.1.rd$rwa1[k] <- mean(dat.f$rwa1[j],na.rm=TRUE)
		dat.f.1.rd$persExtra1[k] <- mean(dat.f$persExtra1[j],na.rm=TRUE)		
		dat.f.1.rd$persAgree1[k] <- mean(dat.f$persAgree1[j],na.rm=TRUE)
		dat.f.1.rd$persConsc1[k] <- mean(dat.f$persConsc1[j],na.rm=TRUE)
		dat.f.1.rd$persNeuro1[k] <- mean(dat.f$persNeuro1[j],na.rm=TRUE)
		dat.f.1.rd$persOpen1[k] <- mean(dat.f$persOpen1[j],na.rm=TRUE)
		dat.f.1.rd$cognition1[k] <- mean(dat.f$cognition1[j],na.rm=TRUE)
		dat.f.1.rd$aggression1[k] <- mean(dat.f$aggression1[j],na.rm=TRUE)
		dat.f.1.rd$age[k] <- mean(dat.f$age[j],na.rm=TRUE)
		dat.f.1.rd$male[k] <- mean(dat.f$male[j],na.rm=TRUE)	
		dat.f.1.rd$christian[k] <- mean(dat.f$christian[j],na.rm=TRUE)
		dat.f.1.rd$white[k] <- mean(dat.f$white[j],na.rm=TRUE)
		dat.f.1.rd$pid1[k] <- mean(dat.f$pid1[j],na.rm=TRUE)
		dat.f.1.rd$ideo1[k] <- mean(dat.f$ideo1[j],na.rm=TRUE)
		dat.f.1.rd$interest1[k] <- mean(dat.f$interest1[j],na.rm=TRUE)
		dat.f.1.rd$educ1[k] <- mean(dat.f$educ1[j],na.rm=TRUE)
		dat.f.1.rd$income1[k] <- mean(dat.f$income1[j],na.rm=TRUE)	
		#Gender composition
		dat.f.1.rd$numF[k] <- length(which(dat.f$male[j]==0))
		dat.f.1.rd$pctF[k] <- (dat.f.1.rd$numF[k]/group.size)
		#Level of group participation (per capita)
		dat.f.1.rd$timesSpoken_rd_g[k] <- sum(dat.f$timesSpoken_rd[j], na.rm=TRUE)
		dat.f.1.rd$timesSpokenPC_rd_g[k] <- dat.f.1.rd$timesSpoken_rd_g[k]/dat.f.1.rd$groupN_rd[k]			
	}
}

#### Descriptive statistics in research design
nrow(dat[which(dat$groupCond=="Individual"),]) #N=760 in individual
nrow(dat[which(dat$groupCond=="Hierarchical"|dat$groupCond=="Horizontal"),]) #N=3213 in group conditions
length(unique(dat$groupID[which(dat$groupCond=="Horizontal")])) #In 406 groups in horizontal
length(unique(dat$groupID[which(dat$groupCond=="Hierarchical")])) #In 365 groups in hierarchical
nrow(dat) #N=3973

## What proportion of group members participated more than once in each deliberation?

#Horizontal
prop.table(table(dat.f$timesSpoken_pt>1)) #76% of group members participate more than once in the chat
prop.table(table(dat.f$timesSpoken_ib>1)) #73% of group members participate more than once in the chat
prop.table(table(dat.f$timesSpoken_rd>1)) #73% of group members participate more than once in the chat

#Hierarchical
#Proportion of leaders who speak at least once in a given module 

prop.table(table(dat.h$timesSpoken_pt[which(dat.h$role_pt=="Leader")] > 1)) #86% of group leaders participate more than once in the chat
prop.table(table(dat.h$timesSpoken_ib[which(dat.h$role_ib=="Leader")] > 1)) #82% of group leaders participate more than once in the chat
prop.table(table(dat.h$timesSpoken_rd[which(dat.h$role_rd=="Leader")] > 1)) #77% of group leaders participate more than once in the chat

# Proportion of advisers who speak at least once in a given module
prop.table(table(dat.h$timesSpoken_pt[which(dat.h$role_pt!="Leader")] > 1)) #80% of group advisers participate more than once in the chat
prop.table(table(dat.h$timesSpoken_ib[which(dat.h$role_ib!="Leader")] > 1)) #76% of group advisers participate more than once in the chat
prop.table(table(dat.h$timesSpoken_rd[which(dat.h$role_rd!="Leader")] > 1)) #74% of group advisers participate more than once in the chat

# Proportion of group members overall in hierarchical condition who speak at least once in a given module
prop.table(table(dat.h$timesSpoken_pt > 1)) #81% of group members participate more than once in the chat
prop.table(table(dat.h$timesSpoken_ib > 1)) #77% of group advisers participate more than once in the chat
prop.table(table(dat.h$timesSpoken_rd > 1)) #74% of group advisers participate more than once in the chat


### Analysis: Prospect theory experiment

## Figure 3: Prospect theory framing effects replicate in groups

B <- 1500
set.seed(43215)

ate.boot <- matrix(NA, nrow=B, ncol=3)
for (i in seq_len(B)){
	resample <- sample(1:nrow(dat.i), nrow(dat.i), replace=TRUE)
	temp <- dat.i[resample,]
	ate.boot[i,1] <- mean(temp$riskyDV_pt[which(temp$lossTrt_pt==1)],na.rm=TRUE) - mean(temp$riskyDV_pt[which(temp$lossTrt_pt==0)],na.rm=TRUE)
	resample <- sample(1:nrow(dat.f.1.pt), nrow(dat.f.1.pt), replace=TRUE)
	temp <- dat.f.1.pt[resample,]
	ate.boot[i,2] <- mean(temp$medDV_pt[which(temp$lossTrt_pt==1)],na.rm=TRUE) - mean(temp$medDV_pt[which(temp$lossTrt_pt==0)],na.rm=TRUE)	
	resample <- sample(1:nrow(dat.h.1), nrow(dat.h.1), replace=TRUE)
	temp <- dat.h.1[resample,]		
	ate.boot[i,3] <- mean(temp$riskyDV_pt[which(temp$lossTrt_pt==1)],na.rm=TRUE) - mean(temp$riskyDV_pt[which(temp$lossTrt_pt==0)],na.rm=TRUE)	
}

ate.df <- data.frame(y=apply(ate.boot,2,mean,na.rm=TRUE), x=as.factor(c("Individual", "Horizontal\ngroup", "Hierarchical\ngroup")), low=apply(ate.boot,2,quantile,0.025,na.rm=TRUE), high=apply(ate.boot,2,quantile,0.975,na.rm=TRUE), low2=apply(ate.boot,2,quantile,0.05,na.rm=TRUE), high2=apply(ate.boot,2,quantile,0.95,na.rm=TRUE))
ate.df$x <- relevel(ate.df$x, ref="Individual")

dev.new(height=5, width=5)
ggplot(ate.df, aes(x,y)) + geom_point() + geom_pointrange(aes(ymin=low, ymax=high)) + geom_linerange(aes(ymin=low2, ymax=high2), lwd=1.5) + geom_hline(yintercept=0, linetype=2) + theme_bw() + xlab("Condition") + ylab("Effect of domain of losses on pr(risky choice)") + ylim(0,0.65) + geom_label(aes(label=round(y,2)), size=2)


## Figure 4: Prospect theory framing effects by horizontal decision rule
set.seed(43215)
ate.boot2 <- matrix(NA, nrow=B, ncol=3)
for (i in seq_len(B)){
	resample <- sample(1:nrow(dat.f.1.pt[which(dat.f.1.pt$un_grp_pt==1),]), nrow(dat.f.1.pt[which(dat.f.1.pt$un_grp_pt==1),]), replace=TRUE)
	temp <- dat.f.1.pt[which(dat.f.1.pt$un_grp_pt==1)[resample],]
	ate.boot2[i,1] <- mean(temp$riskyDV_pt[which(temp$lossTrt_pt==1)],na.rm=TRUE) - mean(temp$riskyDV_pt[which(temp$lossTrt_pt==0)],na.rm=TRUE)
	resample <- sample(1:nrow(dat.f.1.pt), nrow(dat.f.1.pt), replace=TRUE)
	temp <- dat.f.1.pt[resample,]
	ate.boot2[i,2] <- mean(temp$medDV_pt[which(temp$lossTrt_pt==1)],na.rm=TRUE) - mean(temp$medDV_pt[which(temp$lossTrt_pt==0)],na.rm=TRUE)	
	resample <- sample(1:nrow(dat.f.1.pt[which(dat.f.1.pt$maj_grp_pt==1),]), nrow(dat.f.1.pt[which(dat.f.1.pt$maj_grp_pt==1),]), replace=TRUE)
	temp <- dat.f.1.pt[which(dat.f.1.pt$maj_grp_pt==1)[resample],]		
	ate.boot2[i,3] <- mean(temp$riskyDV_pt[which(temp$lossTrt_pt==1)],na.rm=TRUE) - mean(temp$riskyDV_pt[which(temp$lossTrt_pt==0)],na.rm=TRUE)	
}

ate.df2 <- data.frame(y=apply(ate.boot2,2,mean,na.rm=TRUE), x=c("Horizontal\n(Unanimity rule)", "Horizontal\n(Median voter)", "Horizontal\n(Majority rule)"), low=apply(ate.boot2,2,quantile,0.025,na.rm=TRUE), high=apply(ate.boot2,2,quantile,0.975,na.rm=TRUE), low2=apply(ate.boot2,2,quantile,0.05,na.rm=TRUE), high2=apply(ate.boot2,2,quantile,0.95,na.rm=TRUE))

dev.new(height=5, width=5)
ggplot(ate.df2, aes(x,y)) + geom_point() + geom_pointrange(aes(ymin=low, ymax=high)) + geom_linerange(aes(ymin=low2, ymax=high2), lwd=1.5) + geom_hline(yintercept=0, linetype=2) + theme_bw() + xlab("Horizontal decision rule") + ylab("Effect of domain of losses on pr(risky choice)") + ylim(0,0.65) + geom_label(aes(label=round(y,2)), size=2)

## Footnote 49 comparing hierarchical groups to individual decision-makers:

## Comparing hierarchical groups to individual decision-makers, leader-level controls
dat.i.h.pt <- bind_rows(dat.i, dat.h.1.pt)

mod.ih0 <- lm(riskyDV_pt ~ lossTrt_pt * groupCond, dat=dat.i.h.pt)
round(coef(summary(mod.ih0))[nrow(coef(summary(mod.ih0))),4], digits=3) #p < 0.002 without controls
mod.ih1 <- lm(riskyDV_pt ~ lossTrt_pt * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.h.pt)
round(coef(summary(mod.ih1))[nrow(coef(summary(mod.ih1))),4],digits=3) #p < 0.003 with controls at leader-level (demographics, big 5, foreign policy orientations, etc.)

#Comparing hierarchical groups to individual decision-makers, group-level controls
dat.h.1.pt2 <- dat.h.1.pt
a <- grep("mi1", colnames(dat.h.1.pt2))[1] #First remove leader-level covariates so we can merge
b <- grep("smallgrp1", colnames(dat.h.1.pt2))
dat.h.1.pt2 <- dat.h.1.pt2[,-(a:b)]
loc <- grep("_g", colnames(dat.h.1.pt2))
colnames(dat.h.1.pt2)[loc] <- gsub("_g","",colnames(dat.h.1.pt2)[loc])
dat.i.h.pt2 <- bind_rows(dat.i, dat.h.1.pt2)
dat.i.h.pt2$groupCond <- as.factor(dat.i.h.pt2$groupCond)
dat.i.h.pt2$groupCond <- relevel(dat.i.h.pt2$groupCond, ref="Individual")

mod.ih2 <- lm(riskyDV_pt ~ lossTrt_pt * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.h.pt2)
round(coef(summary(mod.ih2))[nrow(coef(summary(mod.ih2))),4],digits=3) #p < 0.002 with controls at group-level

## Statistics in text comparing horizontal groups to individual decision-makers, by decision rule

dat.c.i.f.pt <- bind_rows(dat.i, dat.f.1.pt)
dat.c.i.f.pt$med_DV <- ifelse(is.na(dat.c.i.f.pt$medDV_pt), dat.c.i.f.pt$riskyDV_pt, dat.c.i.f.pt$medDV_pt) #Create common DV for median voter in horizontal, and individuals
dat.c.i.f.pt$un_DV <- ifelse(is.na(dat.c.i.f.pt$un_grp_pt), dat.c.i.f.pt$riskyDV_pt, ifelse(dat.c.i.f.pt$un_grp_pt==1,dat.c.i.f.pt$riskyDV_pt, NA))  #Create common DV for unanimous group in horizontal, and individuals
dat.c.i.f.pt$maj_DV <- ifelse(is.na(dat.c.i.f.pt$maj_grp_pt), dat.c.i.f.pt$riskyDV_pt, ifelse(dat.c.i.f.pt$maj_grp_pt==1,dat.c.i.f.pt$riskyDV_pt, NA))  #Create common DV for majority group in horizontal, and individuals
dat.c.i.f.pt$groupCond <- as.factor(dat.c.i.f.pt$groupCond)
dat.c.i.f.pt$groupCond <- relevel(dat.c.i.f.pt$groupCond, ref="Individual")

#Individual vs horizontal (unanimous groups)
mod.if2 <- lm(un_DV ~ lossTrt_pt * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.c.i.f.pt)
round(coef(summary(mod.if2))[nrow(coef(summary(mod.if2))),4],digits=3) #p < 0.005 with controls 

#Individual vs horizontal (majority-rule groups)
mod.if3 <- lm(maj_DV ~ lossTrt_pt * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.c.i.f.pt)
round(coef(summary(mod.if3))[nrow(coef(summary(mod.if3))),4],digits=3) #p < 0.09 with controls 

# Individual vs horizontal (median voter)
mod.if1 <- lm(med_DV ~ lossTrt_pt * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.c.i.f.pt)
round(coef(summary(mod.if1))[nrow(coef(summary(mod.if1))),4],digits=3) #p < 0.16 with controls

## Figure 5: Intentionality bias replicates in groups
 set.seed(43215)

ate.boot <- matrix(NA, nrow=B, ncol=3)
for (i in seq_len(B)){
	resample <- sample(1:nrow(dat.i), nrow(dat.i), replace=TRUE)
	temp <- dat.i[resample,]
	ate.boot[i,1] <- mean(temp$intentDV_ib[which(temp$fatalTrt_ib==1)],na.rm=TRUE) - mean(temp$intentDV_ib[which(temp$fatalTrt_ib==0)],na.rm=TRUE)
	resample <- sample(1:nrow(dat.f.1.ib), nrow(dat.f.1.ib), replace=TRUE)
	temp <- dat.f.1.ib[resample,]
	ate.boot[i,2] <- mean(temp$medDV_ib[which(temp$fatalTrt_ib==1)],na.rm=TRUE) - mean(temp$medDV_ib[which(temp$fatalTrt_ib==0)],na.rm=TRUE)	
	resample <- sample(1:nrow(dat.h.1), nrow(dat.h.1), replace=TRUE)
	temp <- dat.h.1[resample,]		
	ate.boot[i,3] <- mean(temp$intentDV_ib[which(temp$fatalTrt_ib==1)],na.rm=TRUE) - mean(temp$intentDV_ib[which(temp$fatalTrt_ib==0)],na.rm=TRUE)	
}

ate.df <- data.frame(y=apply(ate.boot,2,mean,na.rm=TRUE), x=as.factor(c("Individual", "Horizontal\ngroup", "Hierarchical\ngroup")), low=apply(ate.boot,2,quantile,0.025,na.rm=TRUE), high=apply(ate.boot,2,quantile,0.975,na.rm=TRUE), low2=apply(ate.boot,2,quantile,0.05,na.rm=TRUE), high2=apply(ate.boot,2,quantile,0.95,na.rm=TRUE))
ate.df$x <- relevel(ate.df$x, ref="Individual")

dev.new(height=5, width=5)
ggplot(ate.df, aes(x,y)) + geom_point() + geom_pointrange(aes(ymin=low, ymax=high)) + geom_linerange(aes(ymin=low2, ymax=high2), lwd=1.5) + geom_hline(yintercept=0, linetype=2) + theme_bw() + xlab("Condition") + ylab("Effect of fatalities on perceived intentionality") + ylim(-0.02,0.25) + geom_label(aes(label=round(y,2)), size=2)

## Figure 6: Intentionality bias effects by horizontal decision rule
set.seed(43215)
ate.boot2 <- matrix(NA, nrow=B, ncol=3)
for (i in seq_len(B)){
	resample <- sample(1:nrow(dat.f.1.ib[which(dat.f.1.ib$un_grp_ib==1),]), nrow(dat.f.1.ib[which(dat.f.1.ib$un_grp_ib==1),]), replace=TRUE)
	temp <- dat.f.1.ib[which(dat.f.1.ib$un_grp_ib==1)[resample],]
	ate.boot2[i,1] <- mean(temp$intentDV_ib[which(temp$fatalTrt_ib==1)],na.rm=TRUE) - mean(temp$intentDV_ib[which(temp$fatalTrt_ib==0)],na.rm=TRUE)
	resample <- sample(1:nrow(dat.f.1.ib), nrow(dat.f.1.ib), replace=TRUE)
	temp <- dat.f.1.ib[resample,]
	ate.boot2[i,2] <- mean(temp$medDV_ib[which(temp$fatalTrt_ib==1)],na.rm=TRUE) - mean(temp$medDV_ib[which(temp$fatalTrt_ib==0)],na.rm=TRUE)	
	resample <- sample(1:nrow(dat.f.1.ib[which(dat.f.1.ib$maj_grp_ib==1),]), nrow(dat.f.1.ib[which(dat.f.1.ib$maj_grp_ib==1),]), replace=TRUE)
	temp <- dat.f.1.ib[which(dat.f.1.ib$maj_grp_ib==1)[resample],]		
	ate.boot2[i,3] <- mean(temp$intentDV_ib[which(temp$fatalTrt_ib==1)],na.rm=TRUE) - mean(temp$intentDV_ib[which(temp$fatalTrt_ib==0)],na.rm=TRUE)	
}

ate.df2 <- data.frame(y=apply(ate.boot2,2,mean,na.rm=TRUE), x=c("Horizontal\n(Unanimity rule)", "Horizontal\n(Median voter)", "Horizontal\n(Majority rule)"), low=apply(ate.boot2,2,quantile,0.025,na.rm=TRUE), high=apply(ate.boot2,2,quantile,0.975,na.rm=TRUE), low2=apply(ate.boot2,2,quantile,0.05,na.rm=TRUE), high2=apply(ate.boot2,2,quantile,0.95,na.rm=TRUE))

dev.new(height=5, width=5)
ggplot(ate.df2, aes(x,y)) + geom_point() + geom_pointrange(aes(ymin=low, ymax=high)) + geom_linerange(aes(ymin=low2, ymax=high2), lwd=1.5) + geom_hline(yintercept=0, linetype=2) + theme_bw() + xlab("Horizontal decision rule") + ylab("Effect of fatalities on perceived intentionality") + ylim(-0.02,0.25) + geom_label(aes(label=round(y,2)), size=2)

## Footnote 50 comparing group conditions to individual decision-makers

# Individual vs horizontal (median voter)
dat.c.i.f.ib <- bind_rows(dat.i, dat.f.1.ib)
dat.c.i.f.ib$med_DV <- ifelse(is.na(dat.c.i.f.ib$medDV_ib), dat.c.i.f.ib$intentDV_ib, dat.c.i.f.ib$medDV_ib) #Create common DV for median voter in horizontal, and individuals
dat.c.i.f.ib$groupCond <- as.factor(dat.c.i.f.ib$groupCond)
dat.c.i.f.ib$groupCond <- relevel(dat.c.i.f.ib$groupCond, ref="Individual")

mod.if0 <- lm(med_DV ~ fatalTrt_ib * groupCond, dat=dat.c.i.f.ib)
round(coef(summary(mod.if0))[nrow(coef(summary(mod.if0))),4],digits=2) #p < 0.94 without controls
mod.if1 <- lm(med_DV ~ fatalTrt_ib * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.c.i.f.ib)
round(coef(summary(mod.if1))[nrow(coef(summary(mod.if1))),4],digits=2) #p < 0.91 with controls

#Individual vs hierarchical: (leader-level controls)
dat.i.h.ib <- bind_rows(dat.i, dat.h.1.ib)
dat.i.h.ib$groupCond <- as.factor(dat.i.h.ib$groupCond)
dat.i.h.ib$groupCond <- relevel(dat.i.h.ib$groupCond, ref="Individual")

mod.ih0 <- lm(intentDV_ib ~ fatalTrt_ib * groupCond, dat=dat.i.h.ib)
round(coef(summary(mod.ih0))[nrow(coef(summary(mod.ih0))),4],digits=2) #p < 0.86 without controls
mod.ih1 <- lm(intentDV_ib ~ fatalTrt_ib * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.h.ib)
round(coef(summary(mod.ih1))[nrow(coef(summary(mod.ih1))),4],digits=2) #p < 0.85 with controls at leader-level

#Individual vs hierarchical: (group-level controls)
dat.h.1.ib2 <- dat.h.1.ib
a <- grep("mi1", colnames(dat.h.1.ib2))[1] #First remove leader-level covariates so we can merge
b <- grep("smallgrp1", colnames(dat.h.1.ib2))
dat.h.1.ib2 <- dat.h.1.ib2[,-(a:b)]
loc <- grep("_g", colnames(dat.h.1.ib2))
colnames(dat.h.1.ib2)[loc] <- gsub("_g","",colnames(dat.h.1.ib2)[loc])
dat.i.h.ib2 <- bind_rows(dat.i, dat.h.1.ib2)
dat.i.h.ib2$groupCond <- as.factor(dat.i.h.ib2$groupCond)
dat.i.h.ib2$groupCond <- relevel(dat.i.h.ib2$groupCond, ref="Individual")

mod.ih2 <- lm(intentDV_ib ~ fatalTrt_ib * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.h.ib2)
round(coef(summary(mod.ih2))[nrow(coef(summary(mod.ih2))),4],digits=2) #p < 0.85 with controls at group-level


### Figure 7: Reactive devaluation experiment displays mixed results
set.seed(43215)

ate.boot <- matrix(NA, nrow=B, ncol=3)
for (i in seq_len(B)){
	resample <- sample(1:nrow(dat.i), nrow(dat.i), replace=TRUE)
	temp <- dat.i[resample,]
	ate.boot[i,1] <- mean(temp$supportDV_rd[which(temp$chinaTrt_rd==1)],na.rm=TRUE) - mean(temp$supportDV_rd[which(temp$chinaTrt_rd==0)],na.rm=TRUE)
	resample <- sample(1:nrow(dat.f.1.rd), nrow(dat.f.1.rd), replace=TRUE)
	temp <- dat.f.1.rd[resample,]
	ate.boot[i,2] <- mean(temp$medDV_rd[which(temp$chinaTrt_rd==1)],na.rm=TRUE) - mean(temp$medDV_rd[which(temp$chinaTrt_rd==0)],na.rm=TRUE)	
	resample <- sample(1:nrow(dat.h.1), nrow(dat.h.1), replace=TRUE)
	temp <- dat.h.1[resample,]		
	ate.boot[i,3] <- mean(temp$supportDV_rd[which(temp$chinaTrt_rd==1)],na.rm=TRUE) - mean(temp$supportDV_rd[which(temp$chinaTrt_rd==0)],na.rm=TRUE)	
}

ate.df <- data.frame(y=apply(ate.boot,2,mean,na.rm=TRUE), x=as.factor(c("Individual", "Horizontal\ngroup", "Hierarchical\ngroup")), low=apply(ate.boot,2,quantile,0.025,na.rm=TRUE), high=apply(ate.boot,2,quantile,0.975,na.rm=TRUE), low2=apply(ate.boot,2,quantile,0.05,na.rm=TRUE), high2=apply(ate.boot,2,quantile,0.95,na.rm=TRUE))
ate.df$x <- relevel(ate.df$x, ref="Individual")

dev.new(height=5, width=5)
ggplot(ate.df, aes(x,y)) + geom_point() + geom_pointrange(aes(ymin=low, ymax=high)) + geom_linerange(aes(ymin=low2, ymax=high2), lwd=1.5) + geom_hline(yintercept=0, linetype=2) + theme_bw() + xlab("Condition") + ylab("Effect of China authorship on support for agreement") + ylim(-0.20,0.25) + geom_label(aes(label=round(y,2)), size=2)

##Figure 8: Reactive devaluation effects by horizontal decision rule
set.seed(43215)
ate.boot2 <- matrix(NA, nrow=B, ncol=3)
for (i in seq_len(B)){
	resample <- sample(1:nrow(dat.f.1.rd[which(dat.f.1.rd$un_grp_rd==1),]), nrow(dat.f.1.rd[which(dat.f.1.rd$un_grp_rd==1),]), replace=TRUE)
	temp <- dat.f.1.rd[which(dat.f.1.rd$un_grp_rd==1)[resample],]
	ate.boot2[i,1] <- mean(temp$supportDV_rd[which(temp$chinaTrt_rd==1)],na.rm=TRUE) - mean(temp$supportDV_rd[which(temp$chinaTrt_rd==0)],na.rm=TRUE)
	resample <- sample(1:nrow(dat.f.1.rd), nrow(dat.f.1.rd), replace=TRUE)
	temp <- dat.f.1.rd[resample,]
	ate.boot2[i,2] <- mean(temp$medDV_rd[which(temp$chinaTrt_rd==1)],na.rm=TRUE) - mean(temp$medDV_rd[which(temp$chinaTrt_rd==0)],na.rm=TRUE)	
	resample <- sample(1:nrow(dat.f.1.rd[which(dat.f.1.rd$maj_grp_rd==1),]), nrow(dat.f.1.rd[which(dat.f.1.rd$maj_grp_rd==1),]), replace=TRUE)
	temp <- dat.f.1.rd[which(dat.f.1.rd$maj_grp_rd==1)[resample],]		
	ate.boot2[i,3] <- mean(temp$supportDV_rd[which(temp$chinaTrt_rd==1)],na.rm=TRUE) - mean(temp$supportDV_rd[which(temp$chinaTrt_rd==0)],na.rm=TRUE)	
}

ate.df2 <- data.frame(y=apply(ate.boot2,2,mean,na.rm=TRUE), x=c("Horizontal\n(Unanimity rule)", "Horizontal\n(Median voter)", "Horizontal\n(Majority rule)"), low=apply(ate.boot2,2,quantile,0.025,na.rm=TRUE), high=apply(ate.boot2,2,quantile,0.975,na.rm=TRUE), low2=apply(ate.boot2,2,quantile,0.05,na.rm=TRUE), high2=apply(ate.boot2,2,quantile,0.95,na.rm=TRUE))

dev.new(height=5, width=5)
ggplot(ate.df2, aes(x,y)) + geom_point() + geom_pointrange(aes(ymin=low, ymax=high)) + geom_linerange(aes(ymin=low2, ymax=high2), lwd=1.5) + geom_hline(yintercept=0, linetype=2) + theme_bw() + xlab("Horizontal decision rule") + ylab("Effect of China authorship on support for agreement") + ylim(-0.20,0.25) + geom_label(aes(label=round(y,2)), size=2)

## Statistics in text comparing individual condition with unanimous horizontal groups

dat.c.i.f.rd <- bind_rows(dat.i, dat.f.1.rd)
dat.c.i.f.rd$un_DV <- ifelse(is.na(dat.c.i.f.rd$un_grp_rd), dat.c.i.f.rd$supportDV_rd, ifelse(dat.c.i.f.rd$un_grp_rd==1,dat.c.i.f.rd$supportDV_rd, NA))  #Create common DV for unanimous group in horizontal, and individuals
dat.c.i.f.rd$groupCond <- as.factor(dat.c.i.f.rd$groupCond)
dat.c.i.f.rd$groupCond <- relevel(dat.c.i.f.rd$groupCond, ref="Individual")

mod.if2 <- lm(un_DV ~ chinaTrt_rd * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.c.i.f.rd)
round(coef(summary(mod.if2))[nrow(coef(summary(mod.if2))),4],digits=2) #p < 0.06 with controls

#### Extensions analysis

### LASSOplus analysis: do certain kinds of leaders help groups avoid these biases?

#Prospect theory experiment
i <- which(is.na(dat.h.1$riskyDV_pt))
y <- dat.h.1$riskyDV_pt[-i]
x <- with(dat.h.1, cbind(mi1, iso1, risk1, govemp1, sdo1, rwa1, persExtra1, persAgree1, persConsc1, persNeuro1, persOpen1, cognition1, aggression1, male, age, educ1, income1, white, pid1, ideo1, interest1))[-i,]
trt <- dat.h.1$lossTrt_pt[-i]

pt.sparse <- sparsereg(y=y, X=x, treat=trt, scale.type="TX", gibbs=1000, burnin=1000, thin=30)
summary(pt.sparse) #The only variable retained is the treatment

#Intentionality bias experiment
i <- which(is.na(dat.h.1$intentDV_ib))
y <- dat.h.1$intentDV_ib[-i]
x <- with(dat.h.1, cbind(mi1, iso1, risk1, govemp1, sdo1, rwa1, persExtra1, persAgree1, persConsc1, persNeuro1, persOpen1, cognition1, aggression1, male, age, educ1, income1, white, pid1, ideo1, interest1))[-i,]
trt <- dat.h.1$fatalTrt_ib[-i]

ib.sparse <- sparsereg(y=y, X=x, treat=trt, scale.type="TX", gibbs=1000, burnin=1000, thin=30)
summary(ib.sparse) #The only variables retained is the treatment

#Reactive devaluation experiment
i <- which(is.na(dat.h.1$supportDV_rd))
y <- dat.h.1$supportDV_rd[-i]
x <- with(dat.h.1, cbind(mi1, iso1, risk1, govemp1, sdo1, rwa1, persExtra1, persAgree1, persConsc1, persNeuro1, persOpen1, cognition1, aggression1, male, age, educ1, income1, white, pid1, ideo1, interest1))[-i,]
trt <- dat.h.1$chinaTrt_rd[-i]

rd.sparse <- sparsereg(y=y, X=x, treat=trt, scale.type="TX", gibbs=1000, burnin=1000, thin=30)
summary(rd.sparse) #No variables retained

### Diversity analysis: Figure 9: More diverse groups are no less susceptible to these three biases

## Prospect theory experiment
f.e <- summary(dat.f.1.pt$diversity_exp)[c(2,5)] #Top quartile: 0.0908
f.d <- summary(dat.f.1.pt$diversity_dem)[c(2,5)] #Top quartile: 0.1348
f.p <- summary(dat.f.1.pt$diversity_disp)[c(2,5)] #Top quartile: 0.06068
f.a <- summary(dat.f.1.pt$diversity_att)[c(2,5)] #Top quartile: 0.0801

mod.h1.pt <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.i)

mod.div.exp.pt <- lm(medDV_pt ~ lossTrt_pt, data=dat.f.1.pt[which(dat.f.1.pt$diversity_exp >= f.e[2]),])
mod.div.dem.pt <- lm(medDV_pt ~ lossTrt_pt, data=dat.f.1.pt[which(dat.f.1.pt$diversity_dem >= f.d[2]),])
mod.div.disp.pt <- lm(medDV_pt ~ lossTrt_pt, data=dat.f.1.pt[which(dat.f.1.pt$diversity_disp >= f.p[2]),])
mod.div.att.pt <- lm(medDV_pt ~ lossTrt_pt, data=dat.f.1.pt[which(dat.f.1.pt$diversity_att >= f.a[2]),])

mod.ldiv.exp.pt <- lm(medDV_pt ~ lossTrt_pt, data=dat.f.1.pt[which(dat.f.1.pt$diversity_exp <= f.e[1]),])
mod.ldiv.dem.pt <- lm(medDV_pt ~ lossTrt_pt, data=dat.f.1.pt[which(dat.f.1.pt$diversity_dem <= f.d[1]),])
mod.ldiv.disp.pt <- lm(medDV_pt ~ lossTrt_pt, data=dat.f.1.pt[which(dat.f.1.pt$diversity_disp <= f.p[1]),])
mod.ldiv.att.pt <- lm(medDV_pt ~ lossTrt_pt, data=dat.f.1.pt[which(dat.f.1.pt$diversity_att <= f.a[1]),])

h.e <- summary(dat.h.1.pt$diversity_exp)[c(2,5)] #Top quartile: 0.08988
h.d <- summary(dat.h.1.pt$diversity_dem)[c(2,5)] #Top quartile: 0.13482
h.p <- summary(dat.h.1.pt$diversity_disp)[c(2,5)] #Top quartile: 0.0615
h.a <- summary(dat.h.1.pt$diversity_att)[c(2,5)] #Top quartile: 0.08101

mod.div.h.exp.pt <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.h.1.pt[which(dat.h.1.pt$diversity_exp >= h.e[2]),])
mod.div.h.dem.pt <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.h.1.pt[which(dat.h.1.pt$diversity_dem >= h.d[2]),])
mod.div.h.disp.pt <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.h.1.pt[which(dat.h.1.pt$diversity_disp >= h.p[2]),])
mod.div.h.att.pt <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.h.1.pt[which(dat.h.1.pt$diversity_att >= h.a[2]),])

mod.ldiv.h.exp.pt <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.h.1.pt[which(dat.h.1.pt$diversity_exp <= h.e[1]),])
mod.ldiv.h.dem.pt <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.h.1.pt[which(dat.h.1.pt$diversity_dem <= h.d[1]),])
mod.ldiv.h.disp.pt <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.h.1.pt[which(dat.h.1.pt$diversity_disp <= h.p[1]),])
mod.ldiv.h.att.pt <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.h.1.pt[which(dat.h.1.pt$diversity_att <= h.a[1]),])

## Intentionality bias experiment
f.e <- summary(dat.f.1.ib$diversity_exp)[c(2,5)] #Top quartile: 0.0887
f.d <- summary(dat.f.1.ib$diversity_dem)[c(2,5)] #Top quartile: 0.1328
f.p <- summary(dat.f.1.ib$diversity_disp)[c(2,5)] #Top quartile: 0.061
f.a <- summary(dat.f.1.ib$diversity_att)[c(2,5)] #Top quartile: 0.0784

mod.h1.ib <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.i)
summary(mod.h1.ib) #0.097 in individual

mod.div.exp.ib <- lm(medDV_ib ~ fatalTrt_ib, data=dat.f.1.ib[which(dat.f.1.ib$diversity_exp >= f.e[2]),])
mod.div.dem.ib <- lm(medDV_ib ~ fatalTrt_ib, data=dat.f.1.ib[which(dat.f.1.ib$diversity_dem >= f.d[2]),])
mod.div.disp.ib <- lm(medDV_ib ~ fatalTrt_ib, data=dat.f.1.ib[which(dat.f.1.ib$diversity_disp >= f.p[2]),])
mod.div.att.ib <- lm(medDV_ib ~ fatalTrt_ib, data=dat.f.1.ib[which(dat.f.1.ib$diversity_att >= f.a[2]),])

mod.ldiv.exp.ib <- lm(medDV_ib ~ fatalTrt_ib, data=dat.f.1.ib[which(dat.f.1.ib$diversity_exp <= f.e[1]),])
mod.ldiv.dem.ib <- lm(medDV_ib ~ fatalTrt_ib, data=dat.f.1.ib[which(dat.f.1.ib$diversity_dem <= f.d[1]),])
mod.ldiv.disp.ib <- lm(medDV_ib ~ fatalTrt_ib, data=dat.f.1.ib[which(dat.f.1.ib$diversity_disp <= f.p[1]),])
mod.ldiv.att.ib <- lm(medDV_ib ~ fatalTrt_ib, data=dat.f.1.ib[which(dat.f.1.ib$diversity_att <= f.a[1]),])

h.e <- summary(dat.h.1.ib$diversity_exp)[c(2,5)] #Top quartile: 0.0899
h.d <- summary(dat.h.1.ib$diversity_dem)[c(2,5)] #Top quartile: 0.1348
h.p <- summary(dat.h.1.ib$diversity_disp)[c(2,5)] #Top quartile: 0.0614
h.a <- summary(dat.h.1.ib$diversity_att)[c(2,5)] #Top quartile: 0.0810

mod.div.h.exp.ib <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.h.1.ib[which(dat.h.1.ib$diversity_exp >= h.e[2]),])
mod.div.h.dem.ib <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.h.1.ib[which(dat.h.1.ib$diversity_dem >= h.d[2]),])
mod.div.h.disp.ib <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.h.1.ib[which(dat.h.1.ib$diversity_disp >= h.p[2]),])
mod.div.h.att.ib <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.h.1.ib[which(dat.h.1.ib$diversity_att >= h.a[2]),])

mod.ldiv.h.exp.ib <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.h.1.ib[which(dat.h.1.ib$diversity_exp <= h.e[1]),])
mod.ldiv.h.dem.ib <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.h.1.ib[which(dat.h.1.ib$diversity_dem <= h.d[1]),])
mod.ldiv.h.disp.ib <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.h.1.ib[which(dat.h.1.ib$diversity_disp <= h.p[1]),])
mod.ldiv.h.att.ib <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.h.1.ib[which(dat.h.1.ib$diversity_att <= h.a[1]),])

## Reactive devaluation experiment
f.e <- summary(dat.f.1.rd$diversity_exp)[c(2,5)] #Top quartile: 0.0861
f.d <- summary(dat.f.1.rd$diversity_dem)[c(2,5)] #Top quartile: 0.1335
f.p <- summary(dat.f.1.rd$diversity_disp)[c(2,5)] #Top quartile: 0.0617
f.a <- summary(dat.f.1.rd$diversity_att)[c(2,5)] #Top quartile: 0.0783

mod.h1.rd <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.i)

mod.div.exp.rd <- lm(medDV_rd ~ chinaTrt_rd, data=dat.f.1.rd[which(dat.f.1.rd$diversity_exp >= f.e[2]),])
mod.div.dem.rd <- lm(medDV_rd ~ chinaTrt_rd, data=dat.f.1.rd[which(dat.f.1.rd$diversity_dem >= f.d[2]),])
mod.div.disp.rd <- lm(medDV_rd ~ chinaTrt_rd, data=dat.f.1.rd[which(dat.f.1.rd$diversity_disp >= f.p[2]),])
mod.div.att.rd <- lm(medDV_rd ~ chinaTrt_rd, data=dat.f.1.rd[which(dat.f.1.rd$diversity_att >= f.a[2]),])

mod.ldiv.exp.rd <- lm(medDV_rd ~ chinaTrt_rd, data=dat.f.1.rd[which(dat.f.1.rd$diversity_exp <= f.e[1]),])
mod.ldiv.dem.rd <- lm(medDV_rd ~ chinaTrt_rd, data=dat.f.1.rd[which(dat.f.1.rd$diversity_dem <= f.d[1]),])
mod.ldiv.disp.rd <- lm(medDV_rd ~ chinaTrt_rd, data=dat.f.1.rd[which(dat.f.1.rd$diversity_disp <= f.p[1]),])
mod.ldiv.att.rd <- lm(medDV_rd ~ chinaTrt_rd, data=dat.f.1.rd[which(dat.f.1.rd$diversity_att <= f.a[1]),])

h.e <- summary(dat.h.1.rd$diversity_exp)[c(2,5)] #Top quartile: 0.077007
h.d <- summary(dat.h.1.rd$diversity_dem)[c(2,5)] #Top quartile: 0.139001
h.p <- summary(dat.h.1.rd$diversity_disp)[c(2,5)] #Top quartile: 0.139001
h.a <- summary(dat.h.1.rd$diversity_att)[c(2,5)] #Top quartile: 0.06427

mod.div.h.exp.rd <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.h.1.rd[which(dat.h.1.rd$diversity_exp >= h.e[2]),])
mod.div.h.dem.rd <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.h.1.rd[which(dat.h.1.rd$diversity_dem >= h.d[2]),])
mod.div.h.disp.rd <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.h.1.rd[which(dat.h.1.rd$diversity_disp >= h.p[2]),])
mod.div.h.att.rd <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.h.1.rd[which(dat.h.1.rd$diversity_att >= h.a[2]),])

mod.ldiv.h.exp.rd <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.h.1.rd[which(dat.h.1.rd$diversity_exp <= h.e[1]),])
mod.ldiv.h.dem.rd <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.h.1.rd[which(dat.h.1.rd$diversity_dem <= h.d[1]),])
mod.ldiv.h.disp.rd <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.h.1.rd[which(dat.h.1.rd$diversity_disp <= h.p[1]),])
mod.ldiv.h.att.rd <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.h.1.rd[which(dat.h.1.rd$diversity_att <= h.a[1]),])


div.out <- data.frame(Type=rep(c("Experience", "Demographics", "Dispositions", "Political Attitudes"),4), Diversity=(rep(c("High", "Low"),each=4)), Rule=rep(c("Horizontal\n(Median voter)", "Hierarchical"),each=8),Module=rep(c("Prospect theory", "Intentionality bias", "Reactive devaluation"), each=16), Y=NA, Low=NA, High=NA, Low2=NA, High2=NA)

div.out[1,5:9] <- ateCalc(mod.div.exp.pt)
div.out[2,5:9] <- ateCalc(mod.div.dem.pt)
div.out[3,5:9] <- ateCalc(mod.div.disp.pt)
div.out[4,5:9] <- ateCalc(mod.div.att.pt)
div.out[5,5:9] <- ateCalc(mod.ldiv.exp.pt)
div.out[6,5:9] <- ateCalc(mod.ldiv.dem.pt)
div.out[7,5:9] <- ateCalc(mod.ldiv.disp.pt)
div.out[8,5:9] <- ateCalc(mod.ldiv.att.pt)
div.out[9,5:9] <- ateCalc(mod.div.h.exp.pt)
div.out[10,5:9] <- ateCalc(mod.div.h.dem.pt)
div.out[11,5:9] <- ateCalc(mod.div.h.disp.pt)
div.out[12,5:9] <- ateCalc(mod.div.h.att.pt)
div.out[13,5:9] <- ateCalc(mod.ldiv.h.exp.pt)
div.out[14,5:9] <- ateCalc(mod.ldiv.h.dem.pt)
div.out[15,5:9] <- ateCalc(mod.ldiv.h.disp.pt)
div.out[16,5:9] <- ateCalc(mod.ldiv.h.att.pt)
div.out[17,5:9] <- ateCalc(mod.div.exp.ib)
div.out[18,5:9] <- ateCalc(mod.div.dem.ib)
div.out[19,5:9] <- ateCalc(mod.div.disp.ib)
div.out[20,5:9] <- ateCalc(mod.div.att.ib)
div.out[21,5:9] <- ateCalc(mod.ldiv.exp.ib)
div.out[22,5:9] <- ateCalc(mod.ldiv.dem.ib)
div.out[23,5:9] <- ateCalc(mod.ldiv.disp.ib)
div.out[24,5:9] <- ateCalc(mod.ldiv.att.ib)
div.out[25,5:9] <- ateCalc(mod.div.h.exp.ib)
div.out[26,5:9] <- ateCalc(mod.div.h.dem.ib)
div.out[27,5:9] <- ateCalc(mod.div.h.disp.ib)
div.out[28,5:9] <- ateCalc(mod.div.h.att.ib)
div.out[29,5:9] <- ateCalc(mod.ldiv.h.exp.ib)
div.out[30,5:9] <- ateCalc(mod.ldiv.h.dem.ib)
div.out[31,5:9] <- ateCalc(mod.ldiv.h.disp.ib)
div.out[32,5:9] <- ateCalc(mod.ldiv.h.att.ib)
div.out[33,5:9] <- ateCalc(mod.div.exp.rd)
div.out[34,5:9] <- ateCalc(mod.div.dem.rd)
div.out[35,5:9] <- ateCalc(mod.div.disp.rd)
div.out[36,5:9] <- ateCalc(mod.div.att.rd)
div.out[37,5:9] <- ateCalc(mod.ldiv.exp.rd)
div.out[38,5:9] <- ateCalc(mod.ldiv.dem.rd)
div.out[39,5:9] <- ateCalc(mod.ldiv.disp.rd)
div.out[40,5:9] <- ateCalc(mod.ldiv.att.rd)
div.out[41,5:9] <- ateCalc(mod.div.h.exp.rd)
div.out[42,5:9] <- ateCalc(mod.div.h.dem.rd)
div.out[43,5:9] <- ateCalc(mod.div.h.disp.rd)
div.out[44,5:9] <- ateCalc(mod.div.h.att.rd)
div.out[45,5:9] <- ateCalc(mod.ldiv.h.exp.rd)
div.out[46,5:9] <- ateCalc(mod.ldiv.h.dem.rd)
div.out[47,5:9] <- ateCalc(mod.ldiv.h.disp.rd)
div.out[48,5:9] <- ateCalc(mod.ldiv.h.att.rd)

ind.out <- data.frame(Type=rep(c("Experience", "Demographics", "Dispositions","Political Attitudes"),each=3), Diversity=NA, Rule=c("Individual"),Module=c("Prospect theory", "Intentionality bias", "Reactive devaluation"), Y=NA, Low=NA, High=NA, Low2=NA, High2=NA)

ind.out[c(1,4,7,10),5:9] <- rep(ateCalc(mod.h1.pt),each=4)
ind.out[c(2,5,8,11),5:9] <- rep(ateCalc(mod.h1.ib),each=4)
ind.out[c(3,6,9,12),5:9] <- rep(ateCalc(mod.h1.rd),each=4)

div.out <- rbind(div.out, ind.out)
div.out$X <- NA
div.out$X[which(div.out$Rule=="Individual")] <- 1
div.out$X[which(div.out$Rule=="Horizontal\n(Median voter)" & div.out$Diversity=="Low")] <- 2
div.out$X[which(div.out$Rule=="Horizontal\n(Median voter)" & div.out$Diversity=="High")] <- 3
div.out$X[which(div.out$Rule=="Hierarchical" & div.out$Diversity=="Low")] <- 4
div.out$X[which(div.out$Rule=="Hierarchical" & div.out$Diversity=="High")] <- 5

div.out$Module <- as.factor(div.out$Module)
div.out$Module <- relevel(div.out$Module, ref="Prospect theory")

dev.new(height=8, width=12)
ggplot(div.out, aes(X, Y)) + geom_point() + geom_pointrange(aes(ymin=Low, ymax=High)) + geom_linerange(aes(ymin=Low2, ymax=High2), lwd=1.5) +  geom_hline(yintercept=0, linetype=2) + theme_bw() + geom_label(aes(label=round(Y,2), fill=Diversity), size=2) + facet_grid(Module ~ Type, scales="free") + scale_x_continuous(breaks=c(1,2.5,4.5), labels=c("Individual","Horizontal\n(Median voter)","Hierarchical")) + xlab("Group condition") + ylab("Average treatment effect")

### Figure 10: More diverse groups are more likely to experience dissension

## Prospect theory experiment

mod.dis.exp.pt <- lm(dissensus ~ lossTrt_pt + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.pt)
mod.dis.dem.pt <- lm(dissensus ~ lossTrt_pt + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.pt)
mod.dis.disp.pt <- lm(dissensus ~ lossTrt_pt + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.pt)
mod.dis.att.pt <- lm(dissensus ~ lossTrt_pt + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.pt)

mod.dis.h.exp.pt <- lm(dissensus ~ lossTrt_pt + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.pt)
mod.dis.h.dem.pt <- lm(dissensus ~ lossTrt_pt + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.pt)
mod.dis.h.disp.pt <- lm(dissensus ~ lossTrt_pt + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.pt)
mod.dis.h.att.pt <- lm(dissensus ~ lossTrt_pt + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.pt)

## Intentionality bias experiment
mod.dis.exp.ib <- lm(dissensus ~ fatalTrt_ib + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.ib)
mod.dis.dem.ib <- lm(dissensus ~ fatalTrt_ib + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.ib)
mod.dis.disp.ib <- lm(dissensus ~ fatalTrt_ib + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.ib)
mod.dis.att.ib <- lm(dissensus ~ fatalTrt_ib + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.ib)

mod.dis.h.exp.ib <- lm(dissensus ~ fatalTrt_ib + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.ib)
mod.dis.h.dem.ib <- lm(dissensus ~ fatalTrt_ib + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.ib)
mod.dis.h.disp.ib <- lm(dissensus ~ fatalTrt_ib + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.ib)
mod.dis.h.att.ib <- lm(dissensus ~ fatalTrt_ib + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.ib)

## Reactive devaluation experiment
mod.dis.exp.rd <- lm(dissensus ~ chinaTrt_rd + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.rd)
mod.dis.dem.rd <- lm(dissensus ~ chinaTrt_rd + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.rd)
mod.dis.disp.rd <- lm(dissensus ~ chinaTrt_rd + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.rd)
mod.dis.att.rd <- lm(dissensus ~ chinaTrt_rd + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.f.1.rd)

mod.dis.h.exp.rd <- lm(dissensus ~ chinaTrt_rd + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.rd)
mod.dis.h.dem.rd <- lm(dissensus ~ chinaTrt_rd + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.rd)
mod.dis.h.disp.rd <- lm(dissensus ~ chinaTrt_rd + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.rd)
mod.dis.h.att.rd <- lm(dissensus ~ chinaTrt_rd + diversity_exp + diversity_dem  + diversity_disp + diversity_att, data=dat.h.1.rd)

dis.out <- data.frame(Type=rep(c("Experience", "Demographics", "Dispositions", "Political Attitudes"),6), Rule=rep(c("Horizontal\n(Median voter)", "Hierarchical"),each=4),Module=rep(c("Prospect theory", "Intentionality bias", "Reactive devaluation"), each=8), Y=NA, Low=NA, High=NA, Low2=NA, High2=NA)

dis.out[1,4:8] <- ateCalc(mod.dis.exp.pt,3)
dis.out[2,4:8] <- ateCalc(mod.dis.dem.pt,4)
dis.out[3,4:8] <- ateCalc(mod.dis.disp.pt,5)
dis.out[4,4:8] <- ateCalc(mod.dis.att.pt,6)
dis.out[5,4:8] <- ateCalc(mod.dis.h.exp.pt,3)
dis.out[6,4:8] <- ateCalc(mod.dis.h.dem.pt,4)
dis.out[7,4:8] <- ateCalc(mod.dis.h.disp.pt,5)
dis.out[8,4:8] <- ateCalc(mod.dis.h.att.pt,6)
dis.out[9,4:8] <- ateCalc(mod.dis.exp.ib,3)
dis.out[10,4:8] <- ateCalc(mod.dis.dem.ib,4)
dis.out[11,4:8] <- ateCalc(mod.dis.disp.ib,5)
dis.out[12,4:8] <- ateCalc(mod.dis.att.ib,5)
dis.out[13,4:8] <- ateCalc(mod.dis.h.exp.ib,3)
dis.out[14,4:8] <- ateCalc(mod.dis.h.dem.ib,4)
dis.out[15,4:8] <- ateCalc(mod.dis.h.disp.ib,5)
dis.out[16,4:8] <- ateCalc(mod.dis.h.att.ib,6)
dis.out[17,4:8] <- ateCalc(mod.dis.exp.rd,3)
dis.out[18,4:8] <- ateCalc(mod.dis.dem.rd,4)
dis.out[19,4:8] <- ateCalc(mod.dis.disp.rd,5)
dis.out[20,4:8] <- ateCalc(mod.dis.att.rd,6)
dis.out[21,4:8] <- ateCalc(mod.dis.h.exp.rd,3)
dis.out[22,4:8] <- ateCalc(mod.dis.h.dem.rd,4)
dis.out[23,4:8] <- ateCalc(mod.dis.h.disp.rd,5)
dis.out[24,4:8] <- ateCalc(mod.dis.h.att.rd,6)

dis.out$X <- NA
dis.out$X[which(dis.out$Rule=="Horizontal\n(Median voter)" & dis.out$Type=="Experience")] <- 3
dis.out$X[which(dis.out$Rule=="Horizontal\n(Median voter)" & dis.out$Type=="Demographics")] <- 1
dis.out$X[which(dis.out$Rule=="Horizontal\n(Median voter)" & dis.out$Type=="Dispositions")] <- 2
dis.out$X[which(dis.out$Rule=="Horizontal\n(Median voter)" & dis.out$Type=="Political Attitudes")] <- 4
dis.out$X[which(dis.out$Rule=="Hierarchical" & dis.out$Type=="Experience")] <- 7
dis.out$X[which(dis.out$Rule=="Hierarchical" & dis.out$Type=="Demographics")] <- 5
dis.out$X[which(dis.out$Rule=="Hierarchical" & dis.out$Type=="Dispositions")] <- 6
dis.out$X[which(dis.out$Rule=="Hierarchical" & dis.out$Type=="Political Attitudes")] <- 8
dis.out$Module <- as.factor(dis.out$Module)
dis.out$Module <- relevel(dis.out$Module, ref="Prospect theory")

dev.new(height=4.5, width=11)
ggplot(dis.out, aes(X, Y)) + geom_point() + geom_pointrange(aes(ymin=Low, ymax=High)) + geom_linerange(aes(ymin=Low2, ymax=High2), lwd=1.5) +  geom_hline(yintercept=0, linetype=2) + theme_bw() + geom_label(aes(label=round(Y,2), fill=Type), size=2) + facet_wrap(~ Module, scales="free") + scale_x_continuous(breaks=c(2,6), labels=c("Horizontal\n(Median voter)","Hierarchical")) + xlab("Group condition") + ylab("Effect of diversity on dissensus")

#### Footnote 73: Effect of gender composition, using a variety of functional forms: 

#PT submodule
#Absolute terms, horizontal groups, linear:

mod.pt.a <- lm(medDV_pt ~ lossTrt_pt * numF, data=dat.f.1.pt)
summary(mod.pt.a) #p < 0.14

#Absolute terms, horizontal groups, factor
mod.pt.a <- lm(medDV_pt ~ lossTrt_pt * as.factor(numF), data=dat.f.1.pt)
summary(mod.pt.a) #p < 0.07 between 0 women groups and 4 women groups, but nothing systematic

#Relative terms, horizontal groups, linear:
mod.pt.r <- lm(medDV_pt ~ lossTrt_pt * pctF, data=dat.f.1.pt)
summary(mod.pt.r) #p < 0.23

#Relative terms, horizontal groups, log:
mod.pt.r <- lm(medDV_pt ~ lossTrt_pt * log(pctF+0.01), data=dat.f.1.pt)
summary(mod.pt.r) #p < 0.23

#IB submodule

#Absolute terms, horizontal groups, linear:

mod.ib.a <- lm(medDV_ib ~ fatalTrt_ib * numF, data=dat.f.1.ib)
summary(mod.ib.a) #p < 0.37

#Absolute terms, horizontal groups, factor
mod.ib.a <- lm(medDV_ib ~ fatalTrt_ib * as.factor(numF), data=dat.f.1.ib)
summary(mod.ib.a) #Nothing systematic

#Relative terms, horizontal groups, linear:
mod.ib.r <- lm(medDV_ib ~ fatalTrt_ib * pctF, data=dat.f.1.ib)
summary(mod.ib.r) #p < 0.68

#Relative terms, horizontal groups, log:
mod.ib.r <- lm(medDV_ib ~ fatalTrt_ib * log(pctF+0.01), data=dat.f.1.ib)
summary(mod.ib.r) #p < 0.66

#RD submodule

#Absolute terms, horizontal groups, linear:

mod.rd.a <- lm(medDV_rd ~ chinaTrt_rd * numF, data=dat.f.1.rd)
summary(mod.rd.a) #p < 0.62

#Absolute terms, horizontal groups, factor
mod.rd.a <- lm(medDV_rd ~ chinaTrt_rd * as.factor(numF), data=dat.f.1.rd)
summary(mod.rd.a) #Nothing systematic

#Relative terms, horizontal groups, linear:
mod.rd.r <- lm(medDV_rd ~ chinaTrt_rd * pctF, data=dat.f.1.rd)
summary(mod.rd.r) #p < 0.97

#Relative terms, horizontal groups, log:
mod.rd.r <- lm(medDV_rd ~ chinaTrt_rd * log(pctF+0.01), data=dat.f.1.rd)
summary(mod.rd.r) #p < 0.43